clear
%给定参数
Nt = 80;
p1 = 100;
p2 = p1/9;
p3 = inf;
h1 = 1000;
h2 = [h1/5, h1/3, h1/2, h1, h1*2, h1*3, h1*5, h1*24];
u = 4*pi*10^(-7);
i = sqrt(-1);
%进行计算
for ix = 1:8
    for  iy = 1: Nt
        w(ix,iy) = 10 ^ ((iy) / 10) / h1;
        lanbta(ix,iy) =  2*pi*sqrt(2) / sqrt(w(ix,iy)*u/p1);
        k1 = sqrt(-i * w(ix,iy) * u / p1);
        k2 = sqrt(-i * w(ix,iy) * u / p2);
        k3 = sqrt(-i * w(ix,iy) * u / p3);
        R2(ix)=tanh((h2(ix) .* k2 + atanh(k2 ./ k3)));
        R1(ix)=tanh( h1 * k1 + atanh((k1 ./ k2) .* R2(ix)));
        P(ix,iy) = R1(ix) .^ 2;
        A(ix,iy) = abs(P(ix,iy));
        B(ix,iy) = angle(R1(ix));
    end

end

figure(9)
loglog(lanbta(1,:)/h1,B(1,:))
hold on
loglog(lanbta(2,:)/h1,B(2,:))
hold on
loglog(lanbta(3,:)/h1,B(3,:))
hold on
loglog(lanbta(4,:)/h1,B(4,:))
hold on
loglog(lanbta(5,:)/h1,B(5,:))
hold on
loglog(lanbta(6,:)/h1,B(6,:))
hold on
loglog(lanbta(7,:)/h1,B(7,:))
hold on
loglog(lanbta(8,:)/h1,B(8,:))
hold off
title("三层大地电磁测深响应曲线 p3 = 100000")
xlabel("lanbta1/h1"); ylabel("相位")
legend("1/5","1/3","1/2","1","2","3","5","24")